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Abstract.  Alexander,  Garcia  and  Alder  [1]  introduced  the  Consistent  Boltzmann  Algorithm  (CBA)  as  a 
simple  variant  of  Direct  Simulation  Monte  Carlo  (DSMC)  for  dense  gases.  In  CBA,  collisions  alter  both 
the  positions  and  velocities  of  the  particles.  In  the  present  work  the  limiting  kinetic  equation  for  CBA  is 
derived;  the  formulation  is  similar  to  the  proof  that  the  Boltzmann  equation  is  the  limiting  kinetic  equation 
for  DSMC  [11].  The  relation  between  this  new  kinetic  equation  and  the  Enskog  equation  is  outlined. 


INTRODUCTION 


Direct  Simulation  Monte  Carlo  (DSMC)  is  presently  the  most  widely  used  numerical  algorithm  in  kinetic 
theory  [3].  The  limiting  kinetic  equation  for  DSMC  is  the  nonlinear  Boltzmann  equation  [11]  so  its  application 
is  restricted  to  dilute  gases.  In  DSMC,  particle  pairs  are  randomly  chosen  to  collide  according  to  the  collision 
probability  for  the  interparticle  potential.  For  example,  for  the  hard  sphere  potential  this  probability  is  propor¬ 
tional  to  the  particles’  relative  speed.  The  post-collision  velocities  are  determined  by  randomly  selecting  the 
collision  angles  and  the  number  of  collisions  on  each  time  step  is  computed  from  the  local  collision  frequency. 
Note  that  in  DSMC  particles  can  be  chosen  to  collide  even  if  their  actual  trajectories  do  not  overlap. 

Recently,  the  Consistent  Boltzmann  Algorithm  (CBA)  was  introduced  as  a  simple  variant  of  DSMC  for  dense 
gases  [1].  Although  CBA  can  be  generalized  to  any  equation  of  state  [2]  here  we  will  only  consider  the  hard 
sphere  gas  with  particle  diameter  a .  In  CBA  the  collision  process  is  as  in  DSMC  with  two  additions.  First, 
when  a  pair  collides  the  unit  apse  vector,  e,  that  is,  the  unit  vector  parallel  to  the  line  connecting  the  centers 
at  impact,  is  computed  from  the  pre-  and  post-collision  velocities  of  the  particles.  Each  particle  is  displaced  a 
distance  a,  one  in  the  direction  e  and  the  other  in  the  direction  — e  (see  Figure  1).  Second,  the  dense  hard 
sphere  collision  frequency,  which  contains  the  so-called  Y -factor,  is  used.  With  these  two  simple  additions  CBA 
yields  the  hard  sphere  equation  of  state  at  all  densities. 

Frezzotti  [5]  and  others  [8]  have  proposed  dense  gas  variants  of  DSMC  based  on  the  Enskog  equation  [9] .  The 
main  advantages  of  CBA  over  Enskog-based  schemes  are  its  simplicity  in  implementation  and  almost  negligible 
effect  on  computational  efficiency  for  a  standard  DSMC  program.  The  transport  coefficients  for  CBA,  obtained 
by  Green-Kubo  analysis,  are  similar  to  those  of  the  Enskog  equation  [1].  As  already  mentioned,  CBA  can  be 
extended  to  potentials  other  than  hard  spheres.  Besides  the  standard  problems  in  kinetic  theory,  CBA  has 
proved  useful  in  the  study  of  granular  materials  [6]  and  nuclear  physics  [7]. 

Until  now  the  principal  deficiency  of  the  Consistent  Boltzmann  Algorithm  was  that  it  lacked  a  complete 
theoretical  foundation.  This  paper  establishes  much  of  that  foundation  by  deriving  the  limiting  kinetic  equation 
for  CBA.  This  equation  is  distinct  from  the  Enskog  equation  and  in  some  respects  is  easier  to  manipulate. 
Besides  its  relation  to  CBA,  this  new  equation  thus  serves  as  a  useful  alternative  to  the  Enskog  equation  in  the 
kinetic  theory  of  dense  gases.  First,  we  give  a  description  of  the  CBA  by  introducing  the  corresponding  Markov 
process.  Second,  we  present  the  limiting  (when  the  number  of  particles  tends  to  infinity)  equation  and  outline 
the  relationship  between  this  equation  and  other  kinetic  equations.  Finally,  for  a  toy  model,  we  solve  the 
limiting  partial  differential  equation  numerically  and  illustrate  the  convergence  behaviour  of  the  corresponding 
stochastic  system. 
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MARKOV  PROCESS 


In  CBA  the  interaction  of  two  particles  {x,  v)  and  {y,  w)  is  determined  by  x*,y*,v*,  and  w*,  the  post-collision 
positions  and  velocities,  respectively.  These  functions  are  defined  as 


v*{x,  v,y,'w,e)  =  v  +  e{e,w  —  v) , 


and 


x*{x,  v,y,w,e)  =  x  +  a 


{v*  —  w*)  —  {v  —  w) 
|(u*  —  w*)  —  {v  —  r(;)||  ’ 


w*{x,  v,y,w,e)  =  w  —  e  (e,  w  —  v) 


y*{x,v,y,w,e)  =  y  -  cr 


{v*  —  w*)  —  {v  —  w) 
||(l!*  —  w*)  —  {v  —  ■ic)||  ’ 


where  e  G  S^,  CTZ^  is  the  unit  sphere,  and  (., .) ,  ||.||  denote  the  scalar  product  and  the  Euclidean  norm  in 
TZ^  ,  respectively.  If  (e,v  —  w)  =0  then  define  x*  =  x  and  y*  =  y .  The  parameter  u  >  0  is  interpreted  as  the 
diameter  of  the  particles.  The  standard  Boltzmann  collision  transformation  is  recovered  in  the  case  cr  =  0 . 

The  Markov  process  related  to  the  collision  step  of  the  CBA  has  states  of  the  form  z  =  {{xi,vi), . . . ,  (a;„,  ?;„)} 
and  the  infinitesimal  generator 


A{^){z)  =  -^  J^^Q{z,i,j,e)^^{J{z,i,j,e))  -  <^{z) 


de . 


(1) 


where  the  jump  transformation  is 


{{xk,Vk)  k^i,j, 
(a:*,u*),if  k  =  i, 
k  =  j, 


and  #  is  a  bounded  measurable  test  function  on  the  state  space.  The  intensity  function  has  the  form 


(2) 


Q{z,i,j,e)  =  Y{p)  h{xi,Xj)B{vi,Vj,e) , 


(3) 


where  is  a  mollifying  kernel  (non-negative  approximation  of  Dirac’s  delta-function)  and  B  is  the  Boltzmann 
collision  kernel.  The  function  h  depends  on  the  cell  structure, 


Hx,  2/)  =  X!  Xc,  (x)  xc,  (y) , 


(4) 


where  \Ci\  denotes  the  Lebesgue  measure  of  the  cell  Ci ,  and  x  is  the  indicator  function.  The  Enskog  Y- 
factor,  Y{p),  gives  the  increase  in  the  collision  frequency  due  to  p,  the  number  density  [9]  in  the  given  cell, 
approximated  by  the  state  2; . 


FIGURE  1.  Two  illustrations  of  CBA  collision  displacement  for  different  apse  vectors.  Before  collision  the  particles 
have  position  and  velocity  (x,v),  (y,w).  After  collision  the  velocities  are  v*  and  w*;  the  shifted  positions  are  indicated 
by  shaded  particles. 
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KINETIC  EQUATION 


The  derivation  of  the  kinetic  equation  for  CBA  roughly  follows  that  presented  in  [11]  for  DSMC.  First  one 
derives  the  equation  satisfied  by  the  limit  (as  n  — >■  oo)  of  the  empirical  measures  of  the  Markov  process.  Then 
the  limiting  equation  for  measure- valued  functions  is  transformed  into  a  form  for  densities.  As  this  derivation 
is  too  lengthy  to  present  here,  we  only  give  the  equation  itself, 


dt 


p{t,x,v)  +  {v,Vx)p{t,x,v)  = 


/  dw  /  deB(v,w,e)  Y(p(t,x*))p{t,x*,v*)p{t,x*,w*) —  Y{p(t,x))p{t,x,v)p(t,x,w)  , 

Jn^  Js^  L 


(5) 


dt‘ 


where  p  is  the  probability  density.  The  mollifier  (4)  is  assumed  to  approach  the  delta-function,  when  the 
number  of  particles  tends  to  infinity. 

Note  that  in  the  case  Y  =  1 ,  a  =  0  the  kinetic  equation  (5)  reduces  to  the  Boltzmann  equation 

^ p{t,x,v)  +  {v,Vx) PitjXjV)  =  f  dw  f  deB{v,w,e)  p{t,x,v*)p{t,x,w*)  —  p{t,x,v)p{t,x,w)  . 
at  j-^s  L  J 

With  5^  =  {e  :  (e,  w  —  v)  >  0}  ,  equation  (5)  takes  the  form 

: p{t,x,v)  +  (VjVx) Pit,x,v)  =  2  /  dw  deB{v,w,e)x  (6) 

•  Jn^  Js\ 

j^y(p(t,  a;  -I-  ae))p{t,x  +  a  e,v*) p(t,x  -t-  ae,w*)  —  Y(p{t,x))  p{t,x,v)  p{t,x,w)^ . 

Equation  (6)  may  be  contrasted  with  the  Enskog  equation  (cf.  [4,  Ch.l6]) 

^p{t,x,v)  +  (v,'Vx)p{t,x,v)  =  2  [  dw  f  deB{v,w,e)x  (7) 

'dt  Jn^  Jsl 

Y  {p{t,  X  +  ^  ere))  p{t,  x,  v*)p{t,  x  ae,  w*)  —  Y  {p{t,  x  —  ^  cr  e))  p{t,  x,  v)  p{t,  x  —  ae,w)  . 

The  revised  Enskog  equation  [10]  is  the  same  as  (7)  with  the  E-factor  replaced  with  the  local-equilibrium 
pair  distribution  function. 


A  TOY  MODEL 


Here  we  study  an  extremely  simplified  model.  For  this  model  the  limiting  partial  differential  equation  is 
solved  numerically,  which  allows  us  to  illustrate  the  convergence  behaviour  of  the  stochastic  system. 

Consider  a  system  where  the  particles  do  not  have  velocities  but  change  their  positions  during  an  interaction. 
The  evolution  of  the  system  is  determined  by  the  generator 


A($)(x) 


2n 


^  /  /i(xi,Xj)[$(J(x,i,j,e))  -  $(x)jde. 


(8) 


where  x  =  (xi, . . . ,  x„)  G  (77^)”  and 


[J(x,i,j,  e)]fc 


X*  ,  if  k^ij, 
Xi  +  a  e  ,  if  k  =  i  , 

Xj  —  a  e ,  if  k  =  j  . 


(9) 


Compared  with  (l)-(3),  the  function  B  is  a  constant,  since  particle  velocities  are  absent.  For  simplicity  we 
take  E  =  1 .  The  analogue  of  the  kinetic  equation  (5)  for  this  toy  model  is  equation 


dt 


p{t,x)  =  B 


x  +  ae)^  —  p{t,x)^ 


de . 


(10) 
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FIGURE  2.  Particle  distribution  in  the  stochastic  systena  (histogram  bars)  and  probability  distribution  from  the 
explicit  difference  scheme  (solid  line)  after  2*®  time  steps. 


FIGURE  3.  Second  and  fourth  moments  versus  time  as  measured  in  the  stochastic  system  (□  and  o)  and  in  the 
explicit  difference  scheme  (x  and  +).  The  solid  and  dashed  lines  go  as  and  respectively. 
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In  the  one-dimensional  case  equation  (10)  takes  the  form 

=  B^p{t,x  +  af  +  p{t,x  -  of  -  2p{t,xfy  (11) 

Note  that  the  unit  sphere  degenerates  to  the  set  {—1;  1} ,  where  each  point  is  given  unit  weight.  With  the 
choice 

D 

B  =  —  ,  for  some  constant  D  >  0 , 
equation  (11)  gives  in  the  limit  ct  — 0  the  partial  differential  equation 

^^p{t,x)=D^p{t,xf .  (12) 

An  alternative  way  of  writing  (12)  is 


d  ,  .  d  ^  d  ,  . 

where  V  =  2Dp{x,t)  can  be  viewed  as  a  nonlinear  diffusion  coefficient. 

The  results  from  numerical  simulations  of  the  stochastic  system  (8),  (9)  and  of  the  explicit  difference  scheme 
for  equation  (11), 


p(t  +  M,x)=p(t.x)  +  AtD  + 


are  shown  in  Figures  2  and  3.  The  stochastic  system  has  5000  particles.  For  both  the  stochastic  system  and 
the  difference  scheme  we  take  a  =  0.3,  D  =  1,  and  At  =  10“^.  The  initial  distribution  is  a  Gaussian  with 
zero  mean  and  unit  variance.  The  distribution  after  a  long  time  is  bullet-shaped,  as  shown  in  Figure  2.  The 
second  and  fourth  moments  go  as  and  as  shown  in  Figure  3.  These  results  may  be  obtained  from 
(12)  using  scaling  arguments.  Note  that  the  distribution  spreads  more  slowly  than  in  the  standard  random 
walk  model  for  which  these  moments  go  as  t  and 
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